1.4 Save dds object

Generate DESeq2 object for differential expression analysis

Here, we will save DESeq2 objects for comparisons within species as well as comparisons across species.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(myTAI)
library(DESeq2)

## Loading required package: S4Vectors
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min
## 
## 
## Attaching package: 'S4Vectors'
## 
## The following objects are masked from 'package:lubridate':
## 
##     second, second<-
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## 
## The following object is masked from 'package:tidyr':
## 
##     expand
## 
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## 
## The following object is masked from 'package:lubridate':
## 
##     %within%
## 
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## The following object is masked from 'package:purrr':
## 
##     reduce
## 
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## 
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## 
## Attaching package: 'MatrixGenerics'
## 
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## 
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## 
## Attaching package: 'Biobase'
## 
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## 
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians

library(IHW)

## 
## Attaching package: 'IHW'
## 
## The following object is masked from 'package:ggplot2':
## 
##     alpha

Within species comparison

Load counts

sample_info_ec <-
  readr::read_csv(file = "data/sample_info_ec.csv")

## Rows: 57 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (9): SampleName, SeqProt, Species, stage_tissue, Stage, Tissue, Experime...
## dbl (1): Replicate
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

sample_info_fd <-
  readr::read_csv(file = "data/sample_info_fd.csv")

## Rows: 34 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (9): SampleName, SeqProt, Species, stage_tissue, Stage, Tissue, Experime...
## dbl (1): Replicate
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

sample_info_fs <-
  readr::read_csv(file = "data/sample_info_fs.csv")

## Rows: 48 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (9): SampleName, SeqProt, Species, stage_tissue, Stage, Tissue, Experime...
## dbl (1): Replicate
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Ec_counts <-
  readr::read_csv("data/Ec_counts.csv")

## Rows: 11582 Columns: 58
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): GeneID
## dbl (57): M_meiospore_low_1, M_meiospore_low_2, M_meiospore_low_3, F_meiospo...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Fd_counts <-
  readr::read_csv("data/Fd_counts.csv")

## Rows: 7907 Columns: 35
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): GeneID
## dbl (34): Fd_f_gametes_1, Fd_f_gametes_2, Fd_f_gametes_3, Fd_m_gametes_1, Fd...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Fs_counts <-
  readr::read_csv("data/Fs_counts.csv")

## Rows: 8291 Columns: 49
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): GeneID
## dbl (48): Fs_F_gametes_1, Fs_F_gametes_2, Fs_F_gametes_3, Fs_M_gametes_1, Fs...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Ec_counts.denoised <-
  readr::read_csv("data/Ec_counts_denoised.csv")

## Rows: 2695 Columns: 58
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): GeneID
## dbl (57): M_meiospore_low_1, M_meiospore_low_2, M_meiospore_low_3, F_meiospo...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Fd_counts.denoised <-
  readr::read_csv("data/Fd_counts_denoised.csv")

## Rows: 2010 Columns: 35
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): GeneID
## dbl (34): Fd_f_gametes_1, Fd_f_gametes_2, Fd_f_gametes_3, Fd_m_gametes_1, Fd...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Fs_counts.denoised <-
  readr::read_csv("data/Fs_counts_denoised.csv")

## Rows: 2214 Columns: 49
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): GeneID
## dbl (48): Fs_F_gametes_1, Fs_F_gametes_2, Fs_F_gametes_3, Fs_M_gametes_1, Fs...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Make into DESeq2 object

Function to reduce the number of lines of code.

# the count data is in tibble format as I have saved it.
create_dds_embyro <- 
  function(count_tibble, metatable, stage_names, rename_list){
    # some error messages
    if(!is_tibble(count_tibble)){return("is not tibble")}
    if(!colnames(count_tibble)[1] == "GeneID"){return("make sure the first column is titled \"GeneID\"")}
    # filter the metatable
    metatable.filt <-
      metatable %>%
      dplyr::filter(Stage %in% stage_names) %>%
      dplyr::arrange(base::match(Stage, stage_names))
    # filter the counts to be loaded as a dds object
    count.embryo <-
      count_tibble %>%
      dplyr::select(
        "GeneID" | 
          metatable.filt$LibName # this would be different for Ectocarpus
      ) %>%
      dplyr::rename_with(~rename_list[.x], -GeneID) %>%
      tibble::remove_rownames() %>%
      tibble::column_to_rownames(var="GeneID")
    # creating a ddsObj for the embryo stages
    simple.model <- 
      as.formula(~Stage)
    ddsObj.raw <- 
      DESeqDataSetFromMatrix(
        round(count.embryo,0),
        colData = metatable.filt,
        design = simple.model)
  }

Fs_embryo_stage_names <- 
  c("24H", "48H", "1w", "3w", "4w")
sample_info_fs[13,2] <- 4
sample_info_fs[14,2] <- 5
sample_info_fs[15,2] <- 6
rename_list <- 
  sample_info_fs %>%
  mutate(new_colname = paste0(Stage, "_", Sex, "_", Replicate)) %>%
  select(LibName, new_colname) %>%
  deframe()

Fs_stage_ddsObj.raw <-
  create_dds_embyro(
    count_tibble = Fs_counts, 
    metatable = sample_info_fs,
    stage_names = Fs_embryo_stage_names,
    rename_list = rename_list)

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

rename_list <- 
  select(sample_info_fd, LibName, Stage, Sex, Replicate) %>%
  mutate(new_colname = paste0(Stage, "_", Sex, "_", Replicate)) %>%
  select(LibName, new_colname) %>%
  deframe()
Fd_embryo_stage_names <- 
  c("E1", "E2", "E3", "E4", "E5", "E6")

Fd_stage_ddsObj.raw <-
  create_dds_embyro(
    count_tibble = Fd_counts, 
    metatable = sample_info_fd,
    stage_names = Fd_embryo_stage_names,
    rename_list = rename_list)

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

# significantly DEG (LRT)
Fs_ddsObj <- 
  DESeq2::DESeq(Fs_stage_ddsObj.raw, test="LRT", reduced=~1)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

Fs_resLRT <- 
  DESeq2::results(Fs_ddsObj)

Fd_ddsObj <- 
  DESeq2::DESeq(Fd_stage_ddsObj.raw, test="LRT", reduced=~1)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

Fd_resLRT <- 
  DESeq2::results(Fd_ddsObj)

# saving with padj filters. With LRT, log foldchange thresholds do not matter.
# see https://hbctraining.github.io/DGE_workshop/lessons/08_DGE_LRT.html
padj_filt = 0.001

tibble::as_tibble(Fs_resLRT, rownames = "GeneID") %>% 
  dplyr::filter(padj < padj_filt)

tibble::as_tibble(Fd_resLRT, rownames = "GeneID") %>% 
  dplyr::filter(padj < padj_filt)

tibble::as_tibble(Fs_resLRT, rownames = "GeneID") %>% 
  ggplot2::ggplot(aes(padj)) +
  ggplot2::geom_histogram() +
  ggplot2::scale_y_log10() +
  ggplot2::labs(
    title = "LRT",
    subtitle = "Fucus serratus",
    x = "padj",
    y = "density (log-scaled)"
  )

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 60 rows containing non-finite values (`stat_bin()`).

tibble::as_tibble(Fd_resLRT, rownames = "GeneID") %>% 
  ggplot2::ggplot(aes(padj)) +
  ggplot2::geom_histogram() +
  ggplot2::scale_y_log10() +
  ggplot2::labs(
    title = "LRT",
    subtitle = "Fucus distichus",
    x = "padj",
    y = "density (log-scaled)"
  )

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 24 rows containing non-finite values (`stat_bin()`).

plot(
  sort(-log10(Fs_resLRT$padj),decreasing = TRUE), 
  xlab="Genes", 
  ylab="-log10(p-value)", 
  main = "Why select 1000 genes?", 
  sub = "Fucus serratus")
abline(v=1000, col="blue")

plot(
  sort(-log10(Fd_resLRT$padj),decreasing = TRUE), 
  xlab="Genes", 
  ylab="-log10(p-value)", 
  main = "Why select 1000 genes?", 
  sub = "Fucus distichus")
abline(v=1000, col="blue")

# select top 1000 DEGs instead.
Fs_LRT.1000 <- 
  tibble::as_tibble(Fs_resLRT, rownames = "GeneID") %>% 
  arrange(padj) %>%
  head(n=1000)
# Max padj:3.828e-40

Fd_LRT.1000 <- 
  tibble::as_tibble(Fd_resLRT, rownames = "GeneID") %>% 
  arrange(padj) %>%
  head(n=1000)
# Max padj:4.405e-37

Save implicated Fucus genes

Fs_LRT.1000 %>%
  readr::write_csv("data/Fs_LRT_1000.csv")
Fd_LRT.1000 %>%
  readr::write_csv("data/Fd_LRT_1000.csv")

6,271 F. serratus and 5,845 F. distichus genes pass the padj filter. This is too much. I have filtered the top 1000 most differentially expressed genes instead by ranking the padj according to the LRT. With the LRT, log2FoldChange is irrelevant.

Furthermore, I have to look a bit more into the LRT. This is because I had thought the LRT tests for genes that can be explained by modelling Stage vs NULL. However, the log2 fold change (MLE): Stage E5 vs E1 makes me worried. However, I think it is fine when ranking with p-values or padj because that is calculated via LRT p-value: '~ Stage' vs '~ 1', which is different to just applying DESeq2::results() on Wald, i.e. Wald test p-value: Stage 4w vs 1w.

Whew.

> Fs_resLRT
log2 fold change (MLE): Stage 4w vs 1w 
LRT p-value: '~ Stage' vs '~ 1' 
DataFrame with 8291 rows and 6 columns

> DESeq2::results(Fs_ddsObj.wald)
log2 fold change (MLE): Stage 4w vs 1w 
Wald test p-value: Stage 4w vs 1w 
DataFrame with 8291 rows and 6 columns

However, I think is fine.

dds for between stages.

This uses the Wald’s test.

summarise_contrasts <- 
  function(
    ddsObj, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,          
    altHypothesis = "greaterAbs",
    stage_names,
    IHW = FALSE,
    output_nsig = TRUE){
    contrast.res <-
      list()
    if(!IHW){
      for (i in 1:(length(stage_names)-1)) {
        contrast.res[[i]] <-
          DESeq2::results(
            ddsObj,
            lfcThreshold = lfcThreshold,
            altHypothesis = altHypothesis,
            contrast = c("Stage", stage_names[i+1], stage_names[i]))
      }
    }
    else if(IHW){
      for (i in 1:(length(stage_names)-1)) {
        contrast.res[[i]] <-
          DESeq2::results(
            ddsObj,
            lfcThreshold = lfcThreshold,
            altHypothesis = altHypothesis,
            filterFun = ihw,
            contrast = c("Stage", stage_names[i+1], stage_names[i]))
      }
    }
    if(!output_nsig){
      return(contrast.res)
    }
    contrast.sig.count <-
      tibble::tibble(
        "contrast" = deseq2_contrasts, 
        "n_sig" = 0, 
        "lfcThreshold" = lfcThreshold,
        "padj_threshold" = padj_threshold,
        "altHypothesis" = altHypothesis,
        "IHW" = IHW
        )
    for (i in 1:length(contrast.res)) {
      contrast.sig.count$n_sig[i] <-
        count(contrast.res[[i]]$padj < padj_threshold, na.rm = T)
    }
    return(contrast.sig.count)
  }

Fs_ddsObj.wald <- 
  DESeq2::DESeq(Fs_stage_ddsObj.raw, test="Wald")

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

Fd_ddsObj.wald <- 
  DESeq2::DESeq(Fd_stage_ddsObj.raw, test="Wald")

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

Fucus serratus

Fs_embryo_stage_names <- 
  c("24H", "48H", "1w", "3w", "4w")
deseq2_contrasts <- c("24H->48H", "48H->1w", "1w->3w", "3w->4w")

Fs_ddsObj.wald_LFC1 <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC2 <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC3 <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC1_IHW <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names,
    IHW = TRUE)
Fs_ddsObj.wald_LFC2_IHW <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names,
    IHW = TRUE)
Fs_ddsObj.wald_LFC3_IHW <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names,
    IHW = TRUE)

Fs_ddsObj.wald_padj0.1 <- 
  base::rbind(Fs_ddsObj.wald_LFC1, 
              Fs_ddsObj.wald_LFC2, 
              Fs_ddsObj.wald_LFC3,
              Fs_ddsObj.wald_LFC1_IHW,
              Fs_ddsObj.wald_LFC2_IHW,
              Fs_ddsObj.wald_LFC3_IHW)

Fs_ddsObj.wald_LFC1 <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC2 <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC3 <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC1_IHW <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names,
    IHW = TRUE)
Fs_ddsObj.wald_LFC2_IHW <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names,
    IHW = TRUE)
Fs_ddsObj.wald_LFC3_IHW <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names,
    IHW = TRUE)

Fs_ddsObj.wald_padj0.05 <- 
  base::rbind(Fs_ddsObj.wald_LFC1, 
              Fs_ddsObj.wald_LFC2, 
              Fs_ddsObj.wald_LFC3,
              Fs_ddsObj.wald_LFC1_IHW,
              Fs_ddsObj.wald_LFC2_IHW,
              Fs_ddsObj.wald_LFC3_IHW)

# significantly nonDEG
Fs_ddsObj.wald_LFC1_less <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC2_less <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC3_less <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)

Fs_ddsObj.wald_less_padj0.1 <- 
  base::rbind(Fs_ddsObj.wald_LFC1_less, 
              Fs_ddsObj.wald_LFC2_less, 
              Fs_ddsObj.wald_LFC3_less)

Fs_ddsObj.wald_LFC1_less <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC2_less <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.wald_LFC3_less <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)

Fs_ddsObj.wald_less_padj0.05 <- 
  base::rbind(Fs_ddsObj.wald_LFC1_less, 
              Fs_ddsObj.wald_LFC2_less, 
              Fs_ddsObj.wald_LFC3_less)

Fucus distichus

Fd_embryo_stage_names <- 
  c("E1", "E2", "E3", "E4", "E5", "E6")

deseq2_contrasts <- c("E1->E2", "E2->E3", "E3->E4", "E4->E5", "E5->E6")

Fd_ddsObj.wald_LFC1 <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC2 <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC3 <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC1_IHW <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names,
    IHW = TRUE)
Fd_ddsObj.wald_LFC2_IHW <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names,
    IHW = TRUE)
Fd_ddsObj.wald_LFC3_IHW <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names,
    IHW = TRUE)

Fd_ddsObj.wald_padj0.1 <- 
  base::rbind(Fd_ddsObj.wald_LFC1, 
              Fd_ddsObj.wald_LFC2, 
              Fd_ddsObj.wald_LFC3,
              Fd_ddsObj.wald_LFC1_IHW,
              Fd_ddsObj.wald_LFC2_IHW,
              Fd_ddsObj.wald_LFC3_IHW)

Fd_ddsObj.wald_LFC1 <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC2 <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC3 <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC1_IHW <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names,
    IHW = TRUE)
Fd_ddsObj.wald_LFC2_IHW <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names,
    IHW = TRUE)
Fd_ddsObj.wald_LFC3_IHW <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names,
    IHW = TRUE)

Fd_ddsObj.wald_padj0.05 <- 
  base::rbind(Fd_ddsObj.wald_LFC1, 
              Fd_ddsObj.wald_LFC2, 
              Fd_ddsObj.wald_LFC3,
              Fd_ddsObj.wald_LFC1_IHW,
              Fd_ddsObj.wald_LFC2_IHW,
              Fd_ddsObj.wald_LFC3_IHW)

# significantly nonDEG
Fd_ddsObj.wald_LFC1_less <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC2_less <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC3_less <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)

Fd_ddsObj.wald_less_padj0.1 <- 
  base::rbind(Fd_ddsObj.wald_LFC1_less, 
              Fd_ddsObj.wald_LFC2_less, 
              Fd_ddsObj.wald_LFC3_less)

Fd_ddsObj.wald_LFC1_less <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC2_less <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.wald_LFC3_less <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)

Fd_ddsObj.wald_less_padj0.05 <- 
  base::rbind(Fd_ddsObj.wald_LFC1_less, 
              Fd_ddsObj.wald_LFC2_less, 
              Fd_ddsObj.wald_LFC3_less)

# combine everything together
Fs_ddsObj.wald_all <-
  base::rbind(Fs_ddsObj.wald_padj0.1,
              Fs_ddsObj.wald_padj0.05,
              Fs_ddsObj.wald_less_padj0.1,
              Fs_ddsObj.wald_less_padj0.05)

Fs_ddsObj.wald_all$contrast <- 
  base::factor(Fs_ddsObj.wald_all$contrast, levels = unique(Fs_ddsObj.wald_all$contrast))

Fd_ddsObj.wald_all <-
  base::rbind(Fd_ddsObj.wald_padj0.1,
              Fd_ddsObj.wald_padj0.05,
              Fd_ddsObj.wald_less_padj0.1,
              Fd_ddsObj.wald_less_padj0.05)

Fd_ddsObj.wald_all$contrast <- 
  base::factor(Fd_ddsObj.wald_all$contrast, levels = unique(Fd_ddsObj.wald_all$contrast))

Plot between-stage significant genes

The difference between performing IHW or not is almost the same.

Fs_ddsObj.wald_all %>%
  dplyr::filter(altHypothesis == "greaterAbs" & IHW == FALSE) %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10() +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus serratus",
    subtitle = "DEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

Fs_ddsObj.wald_all %>%
  dplyr::filter(altHypothesis == "greaterAbs" & IHW == TRUE) %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10() +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus serratus",
    subtitle = "DEG (padj = ihw)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

Fs_ddsObj.wald_all %>%
  dplyr::filter(altHypothesis == "lessAbs") %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10(limits = c(1, 10000)) +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus serratus",
    subtitle = "nonDEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

Fs_ddsObj.wald_all %>%
  dplyr::filter(altHypothesis == "lessAbs") %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::ylim(0, 10000) +
  ggplot2::labs(
    y = "Number of genes",
    x = "Developmental transition",
    title = "Fucus serratus",
    subtitle = "nonDEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

Fd_ddsObj.wald_all %>%
  dplyr::filter(altHypothesis == "greaterAbs" & IHW == FALSE) %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10() +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus distichus",
    subtitle = "DEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

Fd_ddsObj.wald_all %>%
  dplyr::filter(altHypothesis == "greaterAbs" & IHW == TRUE) %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10() +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus distichus",
    subtitle = "DEG (padj = ihw)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

Fd_ddsObj.wald_all %>%
  dplyr::filter(altHypothesis == "lessAbs") %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10(limits = c(1, 10000)) +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus distichus",
    subtitle = "nonDEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

Fd_ddsObj.wald_all %>%
  dplyr::filter(altHypothesis == "lessAbs") %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::ylim(0, 10000) +
  ggplot2::labs(
    y = "Number of genes",
    x = "Developmental transition",
    title = "Fucus distichus",
    subtitle = "nonDEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

denoised dataset (Wald test)

# Fs_embryo_stage_names <- 
#   c("24H", "48H", "1w", "3w", "4w")
# sample_info_fs[13,2] <- 4
# sample_info_fs[14,2] <- 5
# sample_info_fs[15,2] <- 6
rename_list <- 
  sample_info_fs %>%
  mutate(new_colname = paste0(Stage, "_", Sex, "_", Replicate)) %>%
  select(LibName, new_colname) %>%
  deframe()

Fs_stage_ddsObj.denoised.raw <-
  create_dds_embyro(
    count_tibble = Fs_counts.denoised, 
    metatable = sample_info_fs,
    stage_names = Fs_embryo_stage_names,
    rename_list = rename_list)

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

rename_list <- 
  select(sample_info_fd, LibName, Stage, Sex, Replicate) %>%
  mutate(new_colname = paste0(Stage, "_", Sex, "_", Replicate)) %>%
  select(LibName, new_colname) %>%
  deframe()
# Fd_embryo_stage_names <- 
#   c("E1", "E2", "E3", "E4", "E5", "E6")

Fd_stage_ddsObj.denoised.raw <-
  create_dds_embyro(
    count_tibble = Fd_counts.denoised, 
    metatable = sample_info_fd,
    stage_names = Fd_embryo_stage_names,
    rename_list = rename_list)

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

Fs_ddsObj.denoised.wald <- 
  DESeq2::DESeq(Fs_stage_ddsObj.denoised.raw, test="Wald")

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.

## final dispersion estimates

## fitting model and testing

Fd_ddsObj.denoised.wald <- 
  DESeq2::DESeq(Fd_stage_ddsObj.denoised.raw, test="Wald")

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.

## final dispersion estimates

## fitting model and testing

There is no need to use IHW because of the warning Only 1 bin; IHW reduces to Benjamini Hochberg (uniform weights).

Fucus serratus

Fs_embryo_stage_names <- 
  c("24H", "48H", "1w", "3w", "4w")
deseq2_contrasts <- c("24H->48H", "48H->1w", "1w->3w", "3w->4w")

Fs_ddsObj.denoised.wald_LFC1 <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.denoised.wald_LFC2 <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.denoised.wald_LFC3 <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)

Fs_ddsObj.denoised.wald_padj0.1 <- 
  base::rbind(Fs_ddsObj.denoised.wald_LFC1, 
              Fs_ddsObj.denoised.wald_LFC2, 
              Fs_ddsObj.denoised.wald_LFC3)

Fs_ddsObj.denoised.wald_LFC1 <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.denoised.wald_LFC2 <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.denoised.wald_LFC3 <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names)

Fs_ddsObj.denoised.wald_padj0.05 <- 
  base::rbind(Fs_ddsObj.denoised.wald_LFC1, 
              Fs_ddsObj.denoised.wald_LFC2, 
              Fs_ddsObj.denoised.wald_LFC3)

# significantly nonDEG
Fs_ddsObj.denoised.wald_LFC1_less <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.denoised.wald_LFC2_less <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.denoised.wald_LFC3_less <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)

Fs_ddsObj.denoised.wald_less_padj0.1 <- 
  base::rbind(Fs_ddsObj.denoised.wald_LFC1_less, 
              Fs_ddsObj.denoised.wald_LFC2_less, 
              Fs_ddsObj.denoised.wald_LFC3_less)

Fs_ddsObj.denoised.wald_LFC1_less <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.denoised.wald_LFC2_less <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)
Fs_ddsObj.denoised.wald_LFC3_less <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names)

Fs_ddsObj.denoised.wald_less_padj0.05 <- 
  base::rbind(Fs_ddsObj.denoised.wald_LFC1_less, 
              Fs_ddsObj.denoised.wald_LFC2_less, 
              Fs_ddsObj.denoised.wald_LFC3_less)

Fucus distichus

Fd_embryo_stage_names <- 
  c("E1", "E2", "E3", "E4", "E5", "E6")

deseq2_contrasts <- c("E1->E2", "E2->E3", "E3->E4", "E4->E5", "E5->E6")

Fd_ddsObj.denoised.wald_LFC1 <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.denoised.wald_LFC2 <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.denoised.wald_LFC3 <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)

Fd_ddsObj.denoised.wald_padj0.1 <- 
  base::rbind(Fd_ddsObj.denoised.wald_LFC1, 
              Fd_ddsObj.denoised.wald_LFC2, 
              Fd_ddsObj.denoised.wald_LFC3)

Fd_ddsObj.denoised.wald_LFC1 <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.denoised.wald_LFC2 <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.denoised.wald_LFC3 <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names)

Fd_ddsObj.denoised.wald_padj0.05 <- 
  base::rbind(Fd_ddsObj.denoised.wald_LFC1, 
              Fd_ddsObj.denoised.wald_LFC2, 
              Fd_ddsObj.denoised.wald_LFC3)

# significantly nonDEG
Fd_ddsObj.denoised.wald_LFC1_less <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.denoised.wald_LFC2_less <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.denoised.wald_LFC3_less <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.1,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)

Fd_ddsObj.denoised.wald_less_padj0.1 <- 
  base::rbind(Fd_ddsObj.denoised.wald_LFC1_less, 
              Fd_ddsObj.denoised.wald_LFC2_less, 
              Fd_ddsObj.denoised.wald_LFC3_less)

Fd_ddsObj.denoised.wald_LFC1_less <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.denoised.wald_LFC2_less <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 2, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)
Fd_ddsObj.denoised.wald_LFC3_less <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 3, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names)

Fd_ddsObj.denoised.wald_less_padj0.05 <- 
  base::rbind(Fd_ddsObj.denoised.wald_LFC1_less, 
              Fd_ddsObj.denoised.wald_LFC2_less, 
              Fd_ddsObj.denoised.wald_LFC3_less)

# combine everything together
Fs_ddsObj.denoised.wald_all <-
  base::rbind(Fs_ddsObj.denoised.wald_padj0.1,
              Fs_ddsObj.denoised.wald_padj0.05,
              Fs_ddsObj.denoised.wald_less_padj0.1,
              Fs_ddsObj.denoised.wald_less_padj0.05)

Fs_ddsObj.denoised.wald_all$contrast <- 
  base::factor(Fs_ddsObj.denoised.wald_all$contrast, levels = unique(Fs_ddsObj.denoised.wald_all$contrast))

Fd_ddsObj.denoised.wald_all <-
  base::rbind(Fd_ddsObj.denoised.wald_padj0.1,
              Fd_ddsObj.denoised.wald_padj0.05,
              Fd_ddsObj.denoised.wald_less_padj0.1,
              Fd_ddsObj.denoised.wald_less_padj0.05)

Fd_ddsObj.denoised.wald_all$contrast <- 
  base::factor(Fd_ddsObj.denoised.wald_all$contrast, levels = unique(Fd_ddsObj.denoised.wald_all$contrast))

Plot between-stage significant genes

The difference between performing IHW or not is almost the same.

Fs_ddsObj.denoised.wald_all %>%
  dplyr::filter(altHypothesis == "greaterAbs" & IHW == FALSE) %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10() +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus serratus [denoised]",
    subtitle = "DEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

Fs_ddsObj.denoised.wald_all %>%
  dplyr::filter(altHypothesis == "lessAbs") %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10(limits = c(1, 10000)) +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus serratus [denoised]",
    subtitle = "nonDEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

Fs_ddsObj.denoised.wald_all %>%
  dplyr::filter(altHypothesis == "lessAbs") %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::ylim(0, 10000) +
  ggplot2::labs(
    y = "Number of genes",
    x = "Developmental transition",
    title = "Fucus serratus [denoised]",
    subtitle = "nonDEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

Fd_ddsObj.denoised.wald_all %>%
  dplyr::filter(altHypothesis == "greaterAbs" & IHW == FALSE) %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10() +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus distichus [denoised]",
    subtitle = "DEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

Fd_ddsObj.denoised.wald_all %>%
  dplyr::filter(altHypothesis == "lessAbs") %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_y_log10(limits = c(1, 10000)) +
  # ggplot2::annotation_logticks(sides = "l") +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::labs(
    y = "Number of genes (log-scaled)",
    x = "Developmental transition",
    title = "Fucus distichus [denoised]",
    subtitle = "nonDEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

Fd_ddsObj.denoised.wald_all %>%
  dplyr::filter(altHypothesis == "lessAbs") %>%
  ggplot2::ggplot(
    aes(
      x = contrast, 
      y = n_sig, 
      colour = as.factor(lfcThreshold), 
      group = as.factor(lfcThreshold)
      )
    ) +
  ggplot2::geom_line(linewidth = 2, alpha = 0.5) +
  ggplot2::geom_point(size = 3, alpha = 0.5) +
  ggplot2::facet_grid(rows = ggplot2::vars(padj_threshold), labeller = labeller(padj_threshold = label_both)) +
  ggplot2::scale_x_discrete(guide = guide_axis(angle = 30)) +
  ggplot2::ylim(0, 10000) +
  ggplot2::labs(
    y = "Number of genes",
    x = "Developmental transition",
    title = "Fucus distichus [denoised]",
    subtitle = "nonDEG (padj = BH)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position="bottom") +
  ggplot2::guides(colour=guide_legend(title="logFoldChange threshold")) +
  ggsci::scale_color_npg()

Save DEGs

The count isn’t very interesting. So using the same function I made summarise_contrasts() with output_nsig = FALSE. All the data is based on the following thresholds lfcThreshold = 1 and padj_threshold

deseq2_contrasts <- c("24H->48H", "48H->1w", "1w->3w", "3w->4w")
Fs_ddsObj.wald_DEGlist <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names,
    output_nsig = FALSE)
Fs_ddsObj.wald_nonDEGlist <- 
  summarise_contrasts(
    Fs_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names,
    output_nsig = FALSE)

deseq2_contrasts <- c("E1->E2", "E2->E3", "E3->E4", "E4->E5", "E5->E6")
Fd_ddsObj.wald_DEGlist <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names,
    output_nsig = FALSE)
Fd_ddsObj.wald_nonDEGlist <- 
  summarise_contrasts(
    Fd_ddsObj.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names,
    output_nsig = FALSE)

deseq2_contrasts <- c("24H->48H", "48H->1w", "1w->3w", "3w->4w")
Fs_ddsObj.denoised.wald_DEGlist <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fs_embryo_stage_names,
    output_nsig = FALSE)
Fs_ddsObj.denoised.wald_nonDEGlist <- 
  summarise_contrasts(
    Fs_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fs_embryo_stage_names,
    output_nsig = FALSE)

deseq2_contrasts <- c("E1->E2", "E2->E3", "E3->E4", "E4->E5", "E5->E6")
Fd_ddsObj.denoised.wald_DEGlist <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "greaterAbs",
    stage_names = Fd_embryo_stage_names,
    output_nsig = FALSE)
Fd_ddsObj.denoised.wald_nonDEGlist <- 
  summarise_contrasts(
    Fd_ddsObj.denoised.wald, 
    deseq2_contrasts, 
    lfcThreshold = 1, 
    padj_threshold = 0.05,
    altHypothesis = "lessAbs",
    stage_names = Fd_embryo_stage_names,
    output_nsig = FALSE)

Fs_ddsObj.wald_DEGlist %>%
  base::saveRDS("data/Fs_wald_DEG_list.rds")
Fd_ddsObj.wald_DEGlist %>%
  base::saveRDS("data/Fd_wald_DEG_list.rds")

Fs_ddsObj.wald_nonDEGlist %>%
  base::saveRDS("data/Fs_wald_nonDEG_list.rds")
Fd_ddsObj.wald_nonDEGlist %>%
  base::saveRDS("data/Fd_wald_nonDEG_list.rds")

Fs_ddsObj.denoised.wald_DEGlist %>%
  base::saveRDS("data/Fs_denoised_wald_DEG_list.rds")
Fd_ddsObj.denoised.wald_DEGlist %>%
  base::saveRDS("data/Fd_denoised_wald_DEG_list.rds")

Fs_ddsObj.denoised.wald_nonDEGlist %>%
  base::saveRDS("data/Fs_denoised_wald_nonDEG_list.rds")
Fd_ddsObj.denoised.wald_nonDEGlist %>%
  base::saveRDS("data/Fd_denoised_wald_nonDEG_list.rds")

# # to load data
# Fs_ddsObj.wald_DEGlist <-
#   base::readRDS("data/Fs_wald_DEG_list.rds")
# rbind(
#   Fs_ddsObj.wald_DEGlist[[1]] %>%
#     tibble::as_tibble(rownames = "GeneID") %>%
#     dplyr::mutate(contrast = "48H vs 24H"),
#   Fs_ddsObj.wald_DEGlist[[2]] %>%
#     tibble::as_tibble(rownames = "GeneID") %>%
#     dplyr::mutate(contrast = "1w vs 48H"),
#   Fs_ddsObj.wald_DEGlist[[3]] %>%
#     tibble::as_tibble(rownames = "GeneID") %>%
#     dplyr::mutate(contrast = "3w vs 1w"),
#   Fs_ddsObj.wald_DEGlist[[4]] %>%
#     tibble::as_tibble(rownames = "GeneID") %>%
#     dplyr::mutate(contrast = "4w vs 3w"),
#   )

Get session info.

devtools::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.2.2 (2022-10-31)
##  os       macOS Big Sur ... 10.16
##  system   x86_64, darwin17.0
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Europe/Berlin
##  date     2023-11-22
##  pandoc   3.1.6.2 @ /usr/local/bin/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version    date (UTC) lib source
##  annotate               1.76.0     2022-11-01 [1] Bioconductor
##  AnnotationDbi          1.60.2     2023-03-10 [1] Bioconductor
##  Biobase              * 2.58.0     2022-11-01 [1] Bioconductor
##  BiocGenerics         * 0.44.0     2022-11-01 [1] Bioconductor
##  BiocParallel           1.32.6     2023-03-17 [1] Bioconductor
##  Biostrings             2.66.0     2022-11-01 [1] Bioconductor
##  bit                    4.0.5      2022-11-15 [1] CRAN (R 4.2.0)
##  bit64                  4.0.5      2020-08-30 [1] CRAN (R 4.2.0)
##  bitops                 1.0-7      2021-04-24 [1] CRAN (R 4.2.0)
##  blob                   1.2.4      2023-03-17 [1] CRAN (R 4.2.0)
##  bslib                  0.5.1      2023-08-11 [1] CRAN (R 4.2.0)
##  cachem                 1.0.8      2023-05-01 [1] CRAN (R 4.2.0)
##  callr                  3.7.3      2022-11-02 [1] CRAN (R 4.2.0)
##  cli                    3.6.1      2023-03-23 [1] CRAN (R 4.2.0)
##  codetools              0.2-19     2023-02-01 [1] CRAN (R 4.2.0)
##  colorspace             2.1-0      2023-01-23 [1] CRAN (R 4.2.0)
##  crayon                 1.5.2      2022-09-29 [1] CRAN (R 4.2.0)
##  DBI                    1.1.3      2022-06-18 [1] CRAN (R 4.2.0)
##  DelayedArray           0.24.0     2022-11-01 [1] Bioconductor
##  DESeq2               * 1.38.3     2023-01-19 [1] Bioconductor
##  devtools               2.4.5      2022-10-11 [1] CRAN (R 4.2.0)
##  digest                 0.6.33     2023-07-07 [1] CRAN (R 4.2.0)
##  dplyr                * 1.1.3      2023-09-03 [1] CRAN (R 4.2.0)
##  ellipsis               0.3.2      2021-04-29 [1] CRAN (R 4.2.0)
##  evaluate               0.22       2023-09-29 [1] CRAN (R 4.2.2)
##  fansi                  1.0.5      2023-10-08 [1] CRAN (R 4.2.2)
##  farver                 2.1.1      2022-07-06 [1] CRAN (R 4.2.0)
##  fastmap                1.1.1      2023-02-24 [1] CRAN (R 4.2.0)
##  fdrtool                1.2.17     2021-11-13 [1] CRAN (R 4.2.0)
##  forcats              * 1.0.0      2023-01-29 [1] CRAN (R 4.2.0)
##  foreach                1.5.2      2022-02-02 [1] CRAN (R 4.2.0)
##  fs                     1.6.3      2023-07-20 [1] CRAN (R 4.2.0)
##  geneplotter            1.76.0     2022-11-01 [1] Bioconductor
##  generics               0.1.3      2022-07-05 [1] CRAN (R 4.2.0)
##  GenomeInfoDb         * 1.34.9     2023-02-02 [1] Bioconductor
##  GenomeInfoDbData       1.2.9      2023-05-04 [1] Bioconductor
##  GenomicRanges        * 1.50.2     2022-12-16 [1] Bioconductor
##  ggplot2              * 3.4.4      2023-10-12 [1] CRAN (R 4.2.2)
##  ggsci                  3.0.0      2023-03-08 [1] CRAN (R 4.2.0)
##  glue                   1.6.2      2022-02-24 [1] CRAN (R 4.2.0)
##  gtable                 0.3.4      2023-08-21 [1] CRAN (R 4.2.0)
##  hms                    1.1.3      2023-03-21 [1] CRAN (R 4.2.0)
##  htmltools              0.5.6.1    2023-10-06 [1] CRAN (R 4.2.2)
##  htmlwidgets            1.6.2      2023-03-17 [1] CRAN (R 4.2.0)
##  httpuv                 1.6.11     2023-05-11 [1] CRAN (R 4.2.2)
##  httr                   1.4.7      2023-08-15 [1] CRAN (R 4.2.0)
##  IHW                  * 1.26.0     2022-11-01 [1] Bioconductor
##  IRanges              * 2.32.0     2022-11-01 [1] Bioconductor
##  iterators              1.0.14     2022-02-05 [1] CRAN (R 4.2.0)
##  jquerylib              0.1.4      2021-04-26 [1] CRAN (R 4.2.0)
##  jsonlite               1.8.7      2023-06-29 [1] CRAN (R 4.2.0)
##  KEGGREST               1.38.0     2022-11-01 [1] Bioconductor
##  knitr                  1.44       2023-09-11 [1] CRAN (R 4.2.0)
##  labeling               0.4.3      2023-08-29 [1] CRAN (R 4.2.0)
##  later                  1.3.1      2023-05-02 [1] CRAN (R 4.2.2)
##  lattice                0.21-9     2023-10-01 [1] CRAN (R 4.2.2)
##  lifecycle              1.0.3      2022-10-07 [1] CRAN (R 4.2.0)
##  locfit                 1.5-9.8    2023-06-11 [1] CRAN (R 4.2.0)
##  lpsymphony             1.26.3     2023-01-19 [1] Bioconductor (R 4.2.2)
##  lubridate            * 1.9.3      2023-09-27 [1] CRAN (R 4.2.0)
##  magrittr               2.0.3      2022-03-30 [1] CRAN (R 4.2.0)
##  Matrix                 1.5-4.1    2023-05-18 [1] CRAN (R 4.2.0)
##  MatrixGenerics       * 1.10.0     2022-11-01 [1] Bioconductor
##  matrixStats          * 1.0.0      2023-06-02 [1] CRAN (R 4.2.0)
##  memoise                2.0.1      2021-11-26 [1] CRAN (R 4.2.0)
##  mime                   0.12       2021-09-28 [1] CRAN (R 4.2.0)
##  miniUI                 0.1.1.1    2018-05-18 [1] CRAN (R 4.2.0)
##  munsell                0.5.0      2018-06-12 [1] CRAN (R 4.2.0)
##  myTAI                * 1.0.1.9000 2023-10-03 [1] Github (drostlab/myTAI@e159136)
##  pillar                 1.9.0      2023-03-22 [1] CRAN (R 4.2.0)
##  pkgbuild               1.4.2      2023-06-26 [1] CRAN (R 4.2.0)
##  pkgconfig              2.0.3      2019-09-22 [1] CRAN (R 4.2.0)
##  pkgload                1.3.3      2023-09-22 [1] CRAN (R 4.2.0)
##  png                    0.1-8      2022-11-29 [1] CRAN (R 4.2.0)
##  prettyunits            1.2.0      2023-09-24 [1] CRAN (R 4.2.0)
##  processx               3.8.2      2023-06-30 [1] CRAN (R 4.2.0)
##  profvis                0.3.8      2023-05-02 [1] CRAN (R 4.2.0)
##  promises               1.2.1      2023-08-10 [1] CRAN (R 4.2.2)
##  ps                     1.7.5      2023-04-18 [1] CRAN (R 4.2.0)
##  purrr                * 1.0.2      2023-08-10 [1] CRAN (R 4.2.2)
##  R6                     2.5.1      2021-08-19 [1] CRAN (R 4.2.0)
##  RColorBrewer           1.1-3      2022-04-03 [1] CRAN (R 4.2.0)
##  Rcpp                   1.0.11     2023-07-06 [1] CRAN (R 4.2.0)
##  RCurl                  1.98-1.12  2023-03-27 [1] CRAN (R 4.2.0)
##  readr                * 2.1.4      2023-02-10 [1] CRAN (R 4.2.0)
##  remotes                2.4.2.1    2023-07-18 [1] CRAN (R 4.2.2)
##  rlang                  1.1.1      2023-04-28 [1] CRAN (R 4.2.0)
##  rmarkdown              2.25       2023-09-18 [1] CRAN (R 4.2.2)
##  RSQLite                2.3.1      2023-04-03 [1] CRAN (R 4.2.0)
##  rstudioapi             0.15.0     2023-07-07 [1] CRAN (R 4.2.0)
##  S4Vectors            * 0.36.2     2023-02-26 [1] Bioconductor
##  sass                   0.4.7      2023-07-15 [1] CRAN (R 4.2.0)
##  scales                 1.2.1      2022-08-20 [1] CRAN (R 4.2.0)
##  sessioninfo            1.2.2      2021-12-06 [1] CRAN (R 4.2.0)
##  shiny                  1.7.5.1    2023-10-14 [1] CRAN (R 4.2.2)
##  slam                   0.1-50     2022-01-08 [1] CRAN (R 4.2.0)
##  stringi                1.7.12     2023-01-11 [1] CRAN (R 4.2.0)
##  stringr              * 1.5.0      2022-12-02 [1] CRAN (R 4.2.0)
##  SummarizedExperiment * 1.28.0     2022-11-01 [1] Bioconductor
##  tibble               * 3.2.1      2023-03-20 [1] CRAN (R 4.2.0)
##  tidyr                * 1.3.0      2023-01-24 [1] CRAN (R 4.2.0)
##  tidyselect             1.2.0      2022-10-10 [1] CRAN (R 4.2.0)
##  tidyverse            * 2.0.0      2023-02-22 [1] CRAN (R 4.2.0)
##  timechange             0.2.0      2023-01-11 [1] CRAN (R 4.2.0)
##  tzdb                   0.4.0      2023-05-12 [1] CRAN (R 4.2.2)
##  urlchecker             1.0.1      2021-11-30 [1] CRAN (R 4.2.0)
##  usethis                2.2.2      2023-07-06 [1] CRAN (R 4.2.0)
##  utf8                   1.2.3      2023-01-31 [1] CRAN (R 4.2.0)
##  vctrs                  0.6.4      2023-10-12 [1] CRAN (R 4.2.2)
##  vroom                  1.6.4      2023-10-02 [1] CRAN (R 4.2.2)
##  withr                  2.5.1      2023-09-26 [1] CRAN (R 4.2.0)
##  xfun                   0.40       2023-08-09 [1] CRAN (R 4.2.2)
##  XML                    3.99-0.14  2023-03-19 [1] CRAN (R 4.2.0)
##  xtable                 1.8-4      2019-04-21 [1] CRAN (R 4.2.0)
##  XVector                0.38.0     2022-11-01 [1] Bioconductor
##  yaml                   2.3.7      2023-01-23 [1] CRAN (R 4.2.0)
##  zlibbioc               1.44.0     2022-11-01 [1] Bioconductor
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library
## 
## ──────────────────────────────────────────────────────────────────────────────